{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.7 Complex-valued waves\n",
    "====\n",
    "\n",
    "\n",
    "In NGSolve finite element spaces can be built over the complex field and the resulting  complex matrix systems can be solved. This tutorial shows how to compute the solution of the Helmholtz equation with impedance boundary conditions in complex arithmetic. The boundary value problem is to find $u$ satisfying \n",
    "\n",
    "$$\n",
    "-\\Delta u - \\omega^2 u = f\\qquad \\text{ in } \\Omega\n",
    "$$\n",
    "\n",
    "together with the impedance (outgoing) boundary condition\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u }{ \\partial n} - i \\omega u = 0 \n",
    "\\quad \\text{ on } \\partial \\Omega\n",
    "$$\n",
    "\n",
    "where $i =$ `1j` is the imaginary unit.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "from netgen.geom2d import SplineGeometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air = Circle((0.5, 0.5), 0.8).Face()\n",
    "air.edges.name = 'outer'\n",
    "scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()\n",
    "scatterer.edges.name = 'scat'\n",
    "geo = OCCGeometry(air - scatterer, dim=2)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Declare  a complex-valued finite element space "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=5, complex=True)\n",
    "u, v = fes.TnT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = 100\n",
    "pulse = 5e4*exp(-(40**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))\n",
    "Draw(pulse, mesh, order=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forming the system\n",
    "\n",
    "The weak form for $u \\in H^1$:\n",
    "\n",
    "$$\n",
    "\\int_\\Omega\\big[ \\nabla u \\cdot \\nabla \\bar v - \\omega^2 u \\bar v \\big]\n",
    "\\, dx - \n",
    "i \\,\\omega\\, \\int_{\\partial \\Omega} u \\bar v \\, ds = \\int_{\\Omega} f \\bar v\n",
    "$$\n",
    "\n",
    "for all $v$ in $H^1$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forms\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx - omega**2*u*v*dx\n",
    "a += -omega*1j*u*v * ds(\"outer\")\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(pulse * v * dx).Assemble();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, name=\"u\")\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw(gfu, mesh, min=-1, max=1, order=3, animate_complex=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open controls and explore webgui's further viewing options such as viewing real and imaginary parts, viewing absolute value, and viewing with deformation turned on."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
